——— Introducción al análisis espacial ———



Serie de talleres en “Herramientas de análisis cuantitativo y su aplicación en la conservación de la biodiversidad”



Este taller ha sido desarrollado por:

1. Objetivos de aprendizaje

Al final de este taller usted podrá:

  • Identificar las características de los datos espaciales, en sus dos formatos: ráster y vectorial.
  • Conocer los paquetes y funciones en R más utilizados para el manejo de los dos formatos.
  • Utilizar bases de datos de libre acceso disponibles en la web.
  • Integrar y graficar los dos formatos
  • Implemetar algunos análisis espaciales que pueden servir de herramientas para la toma de decisiones en la conservación de la biodiversidad.

2. ¿Qué hacer si no tienes conocimiento previo en R?

Si no tiene concocimientos de R le sugerimos tomar estos cursos que le ayudarán a tener bases y desarrollar sus habilidades.

3. Paquetes en R que vamos a utilizar

  • terra: Análisis espacial de datos (Hijmans et al 2022).

  • raster: Análisis y modelamiento de datos geográficos (Hijmans et al 2022).

  • sf: Simple Features (en inglés) permite el manejo de datos vectoriales (Pebesma et al 2022).

  • stars: permite el manejo de colecciones espacio-temporales y cubos de datos en formato vectorial y raster (Pebesma et al 2022).

  • gdalcubes: permite el manejo de cubos de datos, por ejemplo observaciones de la tierra a través de colecciones de imágenes de satélite (Appel et al 2022).

  • geodata Paquete para decargar datos geogáficos (Hijmans et al 2022).

  • microbenchmark: Mide el tiempo de procesamiento (parece ser mejor que la función system.time() que viene en el paquete base). (microbenchmark)

raster y terra son paquetes en R que permiten el manejo de datos en formato ráster. terra ha sido desarrollado para reemplazar raster, es mucho más rápido y posee mas funcionalidas (entre ellas el manejo de datos vectoriales). En este vínculo puede ver las principales funciones de terra y la comparación de las principales funciones que se encontraban en raster y cómo se llaman ahora en terra.

Por otro lado, el paquete sf (simple features, en inglés) ha reemplazado completamente el paquete sp. Vamos entonces a enfocarnos en sf que permite el manejo de datos vectoriales, como por ejemplo: puntos, líneas y polígonos.

Al final veremos algunas las aplicaciones de los paquetes stars y gdalcubes que permiten el manejo de cubos de datos; como por ejemplo imágenes de satélite que es una colección de varias bandas y años.

Antes de empezar, vamos a instalar (si es necesario) y cargar estos paquetes utilizando la función library().

#install.packages("terra", "raster", "sf", "geodata", "kableExtra", "gdalcubes", "magrittr", "microbenchmark", "stars")
library(terra)
library(raster)
library(sf)
library(geodata)
library(kableExtra)
library(magrittr)
library(gdalcubes)
library(microbenchmark)
library(stars)
remotes::install_github("rlesur/klippy")

4. Para empezar

4.1. Configurar folder de trabajo

Inicie RStudio y abra un nuevo R Script utilizando el menu: File > New File > R Script. Usted podrá copiar los códigos que aquí presentamos y probarlos localmente en su computador. Además habrá unos ejercicios que le permitirán consolidar y aumentar sus habilidades en desarrollar códigos y realizar análisis espaciales.

Digámos que crea un directorio llamado “analisis_espacial” para guardar los archivos de trabajo.

#  "C:/analisis_espacial"

Recuerde que cuando abre RStudio y un nuevo R Script, es importante saber donde esta trabajando con respecto a su disco duro. Para eso puede utilizar la función getwd().

Para asignar el nuevo directorio de trabajo (working directory, en inglés), utilice la función setwd().

# getwd() # consigue el directorio en el cual está trabajando

# setwd("C:/analisis_espacial") 

# en mi caso voy es
  # setwd("C:/Talleres_R/IntroR_espacial")

Ahora que hemos asignado nuestro directorio de trabajo personal, R sabe donde buscar para cargar o guardar archivos!

4.2. Configurar folder de archivos temporales

Cuando se manipulan o se hacen análisis con rásters en R, se generan archivos temporales.

Generalmente uno no sabe esto y menos en donde se guardan estos archivos.

Utilice esta función para ver por defecto en donde se guardan:

#rasterTmpFile()

#Por ejemplo en mi computador es: "C:/Users/jzulo/AppData/Local/Temp/"

Ahora bien, para tener más control de los archivos temnporales es buena practica configurar un folder en su proyecto que se encargue de guardar estos archivos. Al final de su proyecto los podrá borrar. Veremos su utilidad más adelante.

rasterOptions(tmpdir = "./Temp_para_borrar") #Note que estoy utilizando la ruta relativa para crear un folder en mi proyecto "C:/Talleres_R/IntroR_espacial"
rasterTmpFile()
## [1] "C:/Talleres_R/IntroR_espacial/Temp_para_borrar/r_tmp_2024-07-21_194728.352313_28580_92047.grd"

Configurar otros folders para su proyecto

Esto es opcional, pero le permitirá manejar mejor su proyecto.

  • Es bueno crear un folder para guarar todos los datos que bajamos de la red. Por ejemplo creemos un folder llamado “input_data”.
if(!dir.exists("C:/Talleres_R/IntroR_espacial/Input_data")){
dir.create("C:/Talleres_R/IntroR_espacial/Input_data")
}else{
  print("Folder Input_data ya existe")
}
## [1] "Folder Input_data ya existe"
  • Tambien otro folder para guadar sus rasters intermedias que genera dentro de su proyecto pero que probablemente no son los resultados finales.
if(!dir.exists("C:/Talleres_R/IntroR_espacial/Intermediate")){
dir.create("C:/Talleres_R/IntroR_espacial/Intermediate")
}else{
  print("Folder Intermediate ya existe")
}
## [1] "Folder Intermediate ya existe"
  • Finalmente si quiere un folder que contenga los raster finales.
if(!dir.exists("C:/Talleres_R/IntroR_espacial/Outputs")){
dir.create("C:/Talleres_R/IntroR_espacial/Outputs")
}else{
  print("Folder Outputs ya existe")
}
## [1] "Folder Outputs ya existe"

5. ¿Por qué es importante sabe manejar y analizar información espacial?

  • Entender patrones de biodiversidad
  • Hacer prediciones (como porjemplo sobre la distribución de las espcies y el efecto de cambio climático).
  • Planear areas de monitoreo
  • Planeación de la conservación
  • Entre otras

6. ¿Cuáles son los formatos en que encontramos información espacial?

Los principales formatos en el cual encontramos información espacial son: vectorial y ráster.

Los datos vectoriales pueden representar puntos, líneas o polígonos. En general, este formato es común para almacenar y representar observaciones (ej., presencia de especies), objetos lineales (ej., como vías) y áreas geográficas (ej., límites administrativos).

Los datos en ráster (que se concocen también como malla o grilla) almacenan pixeles (celdas). Imagine una fotografía en su cámara digital, o las imágenes de su televisor, o las obtenidas por los satélites, etc. Estas imágenes están compuestas de pixeles.

En adelante utilizaremos el término técnico ráster para hacer referencia a este formato de datos, pero para dar fluidez al texto tambien puede encontrar malla o grilla. Utilizaremos también las palabras pixel o celda para hacer referencia a la unidad mínima de un ráster que contiene información.

7. Ráster

Como mencionamos, el ráster es un formato de datos que almacena imágenes digitales en forma de pixeles (o celdas). En general es una cuadrícula con varias celdas que contienen valores.

Empecemos con unos ejemplos para visualizar y entender mejor de lo que estamos hablando.

7.1. Creando rásters

Vamos a crear un objecto ráster de 10 columnas y 10 filas, utilizando la funcion rast() del paquete terra.

r1 <- rast(ncol=10, 
           nrow=10
           )
r1
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84

Al explorar el objeto r1 vemos lo siguiente:

  • el parámetro class: es SpatRaster (es decir, es básicamente un ráster).
  • el parámetro dimensions: es de 10 filas por 10 columnas (tal y como la definimos) y una capa (nlyr).

Pero además vemos que la función le ha asignado otros atributos como:

  • cood. ref.: la función asigna por defecto coordenadas geográficas lon/lat, usando WGS 84 (World Geodetic System, por sus siglas en inglés).

Ahora bien, la función ha generado únicamente la estructura del ráster, pero todavía el ráster no tiene valores.

Asignemos valores a cada pixel con la funcion values(). Note que la función ncell() calcula el número de pixeles en el ráster. En este caso vamos asignar valores secuenciales.

values(r1) <- 1:ncell(r1)
r1
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :   100

La función values() ha asignado varias cosas:

  • los valores de 1 a 100 (min value: y max value:, respectivamente)
  • extent: por defecto la extensión es en coordenadas geográficas. En este caso longitud (-180 a 180) y latitud (-90 a 90; de polo sur a polo norte).
  • resolution: es la resolución de la grilla, es decir el tamaño de cada pixel o celda, se calcula con base en el número de filas y columnas que indicamos al principio (dimensions:) y con base en el extent:. Es decir, como definimos 10 por 10, la función asigna el siguiente tamaño de pixel 36 x 18 (grados de longitud(x) y latitud(y), respectivamente).
  • el nombre del ráster (name: lyr.1). Lo podemos cambiar más adelante con un nombre que refleje lo que contiene el ráster.
  • el parámetro source: indica la fuente del ráster. En este caso el objeto está en la memoria virtual del computador. Más adelante podrá ver que si descarga un ráster desde su computador podrá utilizar la función sources() para conocer la ruta donde está guardado el ráster.

Visualicemos el ráster que hemos creado, utlizando la función plot() y la función text() para agregar los valores de cada pixel sobre el ráster.

plot(r1) # uilizamos la función plot para visualizar r1
text(r1) # utilizamos la función text para adicionar los valores de cada celda

La grilla que se genera posee valores enteros y la función plot() genera la escala de la barra vertical.

Ya hemos creado nuestro primer ráster!

7.2. Cobertura espacial (extent)

extent: de un ráster representa la cobertura geográfica del ráster, incluyendo celdas sin valor (NAs).

Creemos un ráster basados en el objeto r1 pero reemplacemos los valores de la primera y segunda filas con NA (no values).

r_NA <- r1
r_NA[r_NA < 21]<- NA
r_NA
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :    21 
## max value   :   100

Al visualizar el objeto vemos que las primera dos filas del objeto no se observan (tal y como la creamos), pero la extensión espacial es la misma extent: -180, 180, -90, 90 (xmin, xmax, ymin, ymax)

plot(r_NA, legend=F)
text(r_NA, cex=0.9) # utilizamos la función text para adicionar los valores de cada celda

Comparemos los dos objetos r1 y r_NA.

Como puede ver los dos rásters tiene en mismo extent pero en el segundo no hay datos en las primeras dos filas

Aunque parezca trivial es muy importante tener esto en cuenta, más adelante veremos algunos efectos de esto cuando aplicamos métodos de análisis.

par(mfrow=c(1,2))

# r1
  plot(r1, legend=F) 
  text(r1, cex=0.9) 
  # Adicionemos un punto y las coordenadas en cada esquina de la extensión espacial.
  points(c(xmin(r1), xmax(r1), xmin(r1), xmax(r1)), c(ymax(r1), ymax(r1), ymin(r1), ymin(r1)),               pch=20, col=c("blue", "red"),cex=4)
  # Adicionemos texto
    graphics::text(-40, 97, "xmin=-180; ymax=90", pos=2,  col="blue", cex=1.5)
    graphics::text(180, 97, "xmax=180; ymax=90", pos=2,  col="red", cex=1.5)
    graphics::text(-40, -97, "xmin=-180; ymin=-90", pos=2,  col="blue", cex=1.5)
    graphics::text(180, -97, "xmax=180; ymin=-90", pos=2,  col="red", cex=1.5)

# r_NA
  plot(r_NA, legend=F)
  text(r_NA, cex=0.9) 
 # Adicionemos un punto y las coordenadas en cada esquina de la extensión espacial.
  points(c(xmin(r_NA), xmax(r_NA), xmin(r_NA), xmax(r_NA)), c(ymax(r_NA), ymax(r_NA), ymin(r_NA),            ymin(r1)), pch=20, col=c("blue", "red"),cex=4)
  # Adicionemos texto
    graphics::text(-40, 97, "xmin=-180; ymax=90", pos=2,  col="blue", cex=1.5)
    graphics::text(180, 97, "xmax=180; ymax=90", pos=2,  col="red", cex=1.5)
    graphics::text(-40, -97, "xmin=-180; ymin=-90", pos=2,  col="blue", cex=1.5)
    graphics::text(180, -97, "xmax=180; ymin=-90", pos=2,  col="red", cex=1.5)

7.3. Utilidad de crear un ráster ficticio

El ráster que hemos creado es lo que conocemos como un ráster ficticio (dummy raster en inglés). Las utilidades de construir un ráster ficticio son varias, por ejemplo:

  • podemos utilizar el ráster ficticio para probar funciones/operaciones y verificar los resultados antes de aplicarlas a un ráster grande y con muchos datos. Es decir, podemos verificar que nuestras funciones arrojen los resultados esperados, antes de aplicar dichas funciones a bases de datos reales.

  • podemos compatir en la web nuestros códigos para resolver un problema al cual no le encontramos solución. Es decir compartimos un ejemplo reproducible de lo que queremos hacer y así otras personas pueden sugerir modificaciones. Unas lecturas de la importancia de crear un buen ejemplo reproducible: Elio Campiteli, stackoverplow.com.

Me voy a explicar con un ejemplo para el primer caso.

Primero creemos otro ráster, por ejemplo: r2 será igual a r1+1.

Revise el gráfico de abajo. ¿Es lo esperado?

r2 <- r1+1 # r2 es el resultado de sumarle 1 a r1 creada anteriormente

par(mfrow=c(1,2)) #definamos una ventana de una fila y dos columnas para presentar nuestros rásters

plot(r1, legend=F, main="r1")
text(r1)

plot(r2, legend=F, main="r2")
text(r2)

Creemos ahora otro ráster r3 con base en r1 y r2.

La formula que vamos a utilizar es (r1+r2)/10. Es decir, sumar los valores de r1 y r2 y luego dividirlos entre 10.

En el gráfico de abajo se encuentran los resultados.

¿Es lo esperado?

Parece que sí. Por ejemplo, la primera celda es 0.3, que es el resultado de (1+2)/10 = 0.3

Este es un ejemplo muy simple, pero espero que muestre la utilidad de crear rásters ficticios.

r3 <- (r1+r2)/10

plot(r3, legend=F)
text(r3, digits=2)

7.4. Apilando rásters

Una de las funcionalidades más interesantes con rásters es que las podemos apilar, es decir crear un objeto con un conjunto de capas .

Por ejemplo, apilemos r1, r2, r3 en un objeto que llamaremos r_stack, utilizando la funcion c() del paquete terra (la función análoga en el paquete raster es stack() o brick()).

r_stack <- c(r1,r2,r3)
r_stack
## class       : SpatRaster 
## dimensions  : 10, 10, 3  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## names       : lyr.1, lyr.1, lyr.1 
## min values  :     1,     2,   0.3 
## max values  :   100,   101,  20.1

Vemos que tenemos ahora un objeto con tres capas que se pueden identificar por el nombre (name:en este caso los nombres son iguales) con sus respectivos valores: min values y max values:.

Dentro de este objeto podemos visualizar cada capa de varias formas:

  • La pila de rásters es como una lista de capas, entonces usted puede acceder a cada una utilizando el nombre del objeto seguido de la capa que quiere entre dos parentesis cuadrados, por ejemplo la capa 2 r_stack[[2]].
r_stack[[2]]
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :     2 
## max value   :   101
  • La otra forma es utilizando el signo pesos $ despues del objeto (r_stack$) y le aparecera la lista de las capas que se encuentran en el objeto, identificadas con el nombre (name:). En este caso selecionemos la primera capa (r_stack$lyr.1). El problema es que todas las capas tienen el mismo nombre, entonces usted obtendrá este error:
r_stack$lyr.1
## Error: [subset] you cannot select a layer with a name that is not unique

Aprovechemos y cambiemos los nombres de las capas en esta pila.

names(r_stack) <- c("capa_1", "capa_2", "capa_3") # encadenamos (c) el nombre de cada capa utilizando comillas.

r_stack
## class       : SpatRaster 
## dimensions  : 10, 10, 3  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## names       : capa_1, capa_2, capa_3 
## min values  :      1,      2,    0.3 
## max values  :    100,    101,   20.1

Ahora sí intentemos utilizando el signo $. Parece que si funciona.

r_stack$capa_2
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : capa_2 
## min value   :      2 
## max value   :    101

7.5. Resolución espacial

La resolución espacial de un ráster es básicamente el tamaño del pixel que a su vez representa el área sobre la superficie terrestre.

Por ejemplo en las imágenes Landsat 8 de la (NASA) se muestra la ciudad de Reykjavik en Islandia (July 7, 2019) en donde se pueden apreciar tres resoluciones (30m, 100m y 300m). Es decir desde una resolución fina con más detalle (cada pixel en el terreno es de 30m por 30m) a una más gruesa y borrosa (cada pixel en el terreno es de 300m por 300m).

Como la resolución más fina tiene más detalle, es decir má información, se podrían diferenciar los elementos del terreno y hacer una clasificación detallada y precisa de la superficie terrestre. En la más gruesa es un poco más difícil y se tendría que hacer unas generalizaciones sobre lo que creemos que predomina en cada pixel.

Veamos ahora el concepto de resolución con nuestras rásters ficticias. Creemos otro ráster r4 de 100 columnas y 100 filas y adicionemos el nombre “capa_4”.

r4 <- rast(ncol=100, 
           nrow=100,
           )
values(r4) <- 1:ncell(r4)

names(r4) <- "capa_4" # Note que el objeto es r4  y el nombre de la capa lo cambiamos a "capa_4"

Comparemos r4 con r_stack$capa_1.

r4
## class       : SpatRaster 
## dimensions  : 100, 100, 1  (nrow, ncol, nlyr)
## resolution  : 3.6, 1.8  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : capa_4 
## min value   :      1 
## max value   :  10000
r_stack$capa_1 
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : capa_1 
## min value   :      1 
## max value   :    100

La grilla r4 es mucho mas fina que r_stack$capa_1 porque:

  • la resolución de r4 es de 3.6 por 1.8. Es decir que cada pixel es un rectangulo con un lado de 3.6 y el otro 1.8.

  • y en el ráster r_stack$capa_1 los pixeles son de mayor tamaño (rectangulares tambien) de 36 por 18.

Note que las celdas no necesariamente tienen que ser cuadradas, especialmente cuando tenemos coordenadas geográficas. Más adelante veremos como los rásters projectados sí tienen celdas cuadradas (bueno, en lo posible). Es decir x, y en resolution: deberían tener valores iguales.

Volviendo al ejemplo anterior, podemos también decir que la resolución de r_stack$capa_1 es mucho más gruesa. de hecho 10 veces más gruesa que r4.

Como la resolución de r4 es más fina, entonces hay mas celdas (10,000; dimensions: 100 por 100) que en r_stack$capa_1 (100; dimensions: 10 por 10).

Note que el área seleccionada (extent:) y las coordenadas (coord. ref.:) son las mismas. Visualicemos las dos rásters para ver sus diferencias:

par(mfrow=c(1,2))
plot(r4, legend=F, main="r4, resolución fina")
plot(r_stack$capa_1, legend=F, main="r_stack$capa_1, resolución gruesa")

Ahora bien, ¿Podríamos crear una pila de rásters con área (extent) y coordenadas (coord. ref.) similares, pero resoluciones y dimensiones diferentes?

Veamos:

rstack_new <- c(r_stack$capa_1, r4)
## Error: [rast] number of rows and/or columns do not match

Pues NO, el Error nos dice que el número de columnas y filas (dimensions) no coinciden.

7.6. ¿Qué más podemos hacer con rásters?

Ya vimos algunas cosas que podemos hacer con rásters, como por ejemplo apilarlas. Pero hay algunos métodos importantes cuando se quieren hacer análisis utilizando rásters. Cuando nos referimos a métodos hacemos referencia a una serie de funciones que se pueden aplicar a uno o varios rásters para obtener información relevante en su investigación.

No podemos en este taller ver todas las funciones, pero al final usted tendrá la capacidad de explorar y aplicar más funciones.

Por ejemplo, ya vimos cómo podemos crear un ráster utilizando la función rast(), apilar varios rásters con c(), y cómo seleccionar algunas capas dentro de una pila de rásters utilizando [[]] o $.

Ahora veamos algunos métodos y funciones que puden ser mucha utilidad en el análisis espacial.

7.6.1. Métodos locales

Los métodos locales hacen referencia a hacer cálculos entre celdas que se sobreponen en una pila de rásters. Es decir la localización espacial de las celdas es la misma. La figura muestra cómo podríamos calcular la media de dos rásters que se sobreponen.

Hay varias funciones para hacer cálculos de este tipo. Por ahora utilizaremos la función app() para demostrar este método.

Veamos el código en R para calcular la media de la pila de tres capas que hemos creado; es decir la funcion app() en terra toma los tres valores de cada pixel de cada capa que se sobreponen y calcula la media (mean). Esto se hace en toda la extensión del área que definimos.

r_stack_media <- app(x = r_stack, # Este es el objeto que representa la pila de rásters
                     fun = mean   # Esta es la función para calcular la media 
                     )
r_stack_media
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : mean 
## min value   :  1.1 
## max value   : 73.7

Note que el nuevo objeto r_stack_media representa la media de las tres capas en r_stack. R asigna el nombre de la función a la capa final en esta caso name: es mean.

Visualicemos el promedio de esto tres rásters.

plot(r_stack_media)
text(r_stack_media, digits=1, cex=0.8)

Tips!

Note que los decimales están indicados con un punto ., para los miles se utiliza la coma ,. Los países del norte utilizan esta nomenclatura y R tambien.

En Colombia es al contrario (no sabría decir si en toda latinoamérica es igual). Por qué sucede esto? No lo sé.

Lo que sé es que usted tendrá que ajustar las bases de datos para que tengan el mismo formato requerido en R. En especial tablas en excel que le gustaría incorporarlas a un análisis espacial.

Usted ya puede ver las potenciales aplicaciones de estas funciones en un pila de rásters. Por ejemplo, calcular la media o la varianza de la temperatura durante un año (una pila de 12 capas) en su área de estudio.

Ahora veamos un ejemplo utilizando datos globales reales!!!!

Datos globales y algunas aplicaciones

Vamos ahora a manejar algunos bases de datos de acceso gratuito (open access, en inglés) que pueden ser utilizadas en análisis espaciales.

Por ejemplo, los datos globales de clima, como temperatura y precipitación, pueden ser utilizados ampliamente para entender los efectos del cambio climático en la distribución de las especies.

WorldClim y CHELSA son dos iniciativas que ofrecen datos de buena calidad a nivel global.

Existen además algunos paquetes en R que permiten descargar esta bases de datos. Por ejemplo, en el paquete raster existe la función getData() que permite descargar datos de WorldClim. Sin embargo, esta función será reemplazada por el paquete geoData. Veamos un ejemplo

Temperatura media anual a nivel global (WorldClim)

Descarguemos la temperatura media anual a nivel global a una resolución gruesa (10 minutos que es ~340 km2) desde WorldClim. Tenga en cuenta que descargará un archivo comprimido (.zip) a su computador.

Existe tambien la posibilidad de descargar las 19 bioclimatic variables disponibles utilizando var = ‘bio’. Por ahora descargemos la temperatura media (tagv):

tmean <- geodata::worldclim_global(
    var="tavg", # temperatura media
    res=10,     # resolución gruesa
    path = "C:/Talleres_R/IntroR_espacial" # folder donde desea guardar los datos
                      )

Puede ver que son 12 capas (una por cada mes, ver el atributo nlyr en dimensions:).

tmean
## class       : SpatRaster 
## dimensions  : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : wc2.1_10m_tavg_01.tif  
##               wc2.1_10m_tavg_02.tif  
##               wc2.1_10m_tavg_03.tif  
##               ... and 9 more source(s)
## names       : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ... 
## min values  :    -45.8840,   -44.80000,   -57.92575,   -64.19250,   -64.81150,   -64.35825, ... 
## max values  :     34.0095,    32.82425,    32.90950,    34.19375,    36.25325,    38.35550, ...

Visualicemos dos capas (Enero y Agosto)

par(mfrow=c(1,2)) 

terra::plot(tmean$wc2.1_10m_tavg_01, main="Temperatura media - Enero", range=c(-50, 40), mar=c(4, 3, 4, 3), plg=list(shrink=0.9, cex=.8), 
    pax=list(side=1:2, cex.axis=.6))
terra::plot(tmean$wc2.1_10m_tavg_10, main="Temperatura media - Agosto", range=c(-50, 40), mar=c(4, 2, 4, 4), legend=FALSE, 
    pax=list(side=c(1,4), cex.axis=.6))

Excelente, ya vemos algunos rásters con datos reales. Ahora tratemos de aplicar algunas funciones.

Funciones locales

Vamos a aplicar las funciones stdev()y app() al objeto tmean que tiene 12 capas.

  • Medidas de dispersion stdev()

Veamos en donde hay mayor variación de la temperatura durante el año, utilizando la función stdv() del paquete terra.

tmean_sd <- stdev(x = tmean, # SpatRaster object
                  pop = T, # Si es verdadero (T), entonces la desviación estándard de la población es calculada
                  na.rm = T # Si es verdadero (T), entonces valores NA son ignorados
                  )

¿Es lo esperado?

plot(tmean_sd)

Tips!

Recuerde que en RStudio puede iluminar la función (ej., stdev) y oprimir F1 para obtener ayuda, es decir información sobre los argumentos que se encuentran dentro de una función.

  • Calcular la temperatura media utilizando app()

También podemos calcular la temperatura media del año, utilizando la función app() del paquete terra.

tmean_mean <- terra::app(x = tmean,  # SpatRaster
                  fun = mean) # función
plot(tmean_mean, main="Temperatura media anual (°C)")

7.6.2. Métodos zonales y globales

Los métodos zonales y globales hacen referencia a hacer cálculos del área total del ráster utilizando la función global() o zonas de un ráster utilizando la función zonal() (estas zonas se pueden definir con otro capa). Por ejemplo, la figura muestra cómo podemos calcular la sumatoria de una zona en el ráster usando un polígono.

Veamos un ejemplo de cómo podemos utilizar la base de datos global de la huella humana (Human footprint, en inglés) para calcular el promedio de huella humana para cada país del mundo.

La huella humana reune la mayor parte de las acciones humanas con el potencial de deteriorar la naturaleza (ver el trabajo de Venter et al 2016 para mayores detalles).

La escala es 0 = muy baja, hasta 50 = muy alta huella humana.

Note que en nuestros códigos vamos a utilizar el simbolo %>%, conocido como forward-pipe operator en inglés, que envia el resultado de una expresión a la otra función dentro de un mismo objeto. Este operador%>% fue originalmente desarrollado paquete magrittr.

En este caso cargamos la base de datos “footprint” y después la agregamos 10 veces utilizando la función aggregate(). Tenga en cuenta que la función aggregate() se encuentra en varios paquetes, entonces es buena práctica especificar que estámos utilizando la función del paquete terra, asi: terra:aggregate().

¿Por qué agregamos el ráster huella humana? Porque queremos la huella humana en un ráster de resolución más gruesa para acelerar el procesamiento de nuestros siguientes pasos (es decir menos tiempo de procesamiento). Más adelante veremos unos ejemplos del tiempo de procesamiento en diferentes resoluciones.

huella <- geodata::footprint(year = 2009,
                             path ="C:/Talleres_R/IntroR_espacial")%>% 
          terra::aggregate(fact=10, fun=mean)
## 
|---------|---------|---------|---------|
=========================================
                                          
huella
## class       : SpatRaster 
## dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : wildareas-v3-2009-human-footprint 
## min value   :                              0.00 
## max value   :                             45.21

Ahora carguemos los países del mundo. Es una base de datos que contiene información vectorial. Más adelante dedicaremos una sección a este formato.

paises <- geodata::world(resolution = 5,
                         level = 0,
                        path = "C:/Talleres_R/IntroR_espacial")
paises
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 231, 2  (geometries, attributes)
##  extent      : -180, 180, -90, 83.65625  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : GID_0      NAME_0
##  type        : <chr>       <chr>
##  values      :   ABW       Aruba
##                  AFG Afghanistan
##                  AGO      Angola

Ahora si, apliquemos la función zonal() del paquete terra para calcular el promedio de la huella humana en cada país del mundo.

huella_paises <- terra::zonal(huella,
                              paises,
                              as.raster=T, # si es F obtiene solo los datos.
                              fun = mean
                              )
plot(huella_paises)

7.6.3. Seleccionado un área geográfica dentro de un ráster

Es muy común que usted tenga un ráster con una extensión global, pero usted quiera enfocar su análisis en una región específica. Para esto vamos a utilizar la función crop() del paquete terra::.

Por ejemplo, seleccionemos los datos de la huella humana solo para Colombia.

Por ahora transformemos la información vectorial a ráster, utilizando la función rasterize() del paquete terra::.

Note que vamos a utilizar el ráster huella como plantilla para rasterize() los polígonos paises; es decir obtendremos un ráster de las mismas dimensions:, resolution:, extent y coord. ref.:.

Para los valores del ráster utilizaremos el nombre de la columna en el objeto vectorial que contiene el nombre de los países, es decir NAME_0.

paises_r <- terra::rasterize(paises, huella, field="NAME_0")
paises_r
## class       : SpatRaster 
## dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## categories  : NAME_0 
## name        :      NAME_0 
## min value   : Afghanistan 
## max value   :       Åland
plot(paises_r)

Las categorias del objeto paises_r pueden ser inspeccionadas con la función levels() del paquete terra. Generalmente son dos columnas, la primera es un numero entero que identifica los valores de la celda y la segunda es una etiqueta que indica el nombre de la categoria. Si revisa la tabla podrá ver que Colombia es el valor 43.

knitr::kable(head(terra::levels(paises_r)))%>%
  kableExtra::kable_paper("hover", full_width = F) %>%
  kableExtra::scroll_box(width = "100%", height = "200px")
ID NAME_0
0 Afghanistan
1 Akrotiri and Dhekelia
2 Albania
3 Algeria
4 American Samoa
5 Andorra
6 Angola
7 Antarctica
8 Antigua and Barbuda
9 Argentina
10 Armenia
11 Aruba
12 Australia
13 Austria
14 Azerbaijan
15 Bahamas
16 Bahrain
17 Bangladesh
18 Barbados
19 Belarus
20 Belgium
21 Belize
22 Benin
23 Bhutan
24 Bolivia
25 Bonaire, Saint Eustatius and Saba
26 Bosnia and Herzegovina
27 Botswana
28 Brazil
29 Brunei
30 Bulgaria
31 Burkina Faso
32 Burundi
33 Cabo Verde
34 Cambodia
35 Cameroon
36 Canada
37 Caspian Sea
38 Cayman Islands
39 Central African Republic
40 Chad
41 Chile
42 China
43 Colombia
44 Comoros
45 Congo
46 Cook Islands
47 Costa Rica
48 Croatia
49 Cuba
50 Curaçao
51 Cyprus
52 Czech Republic
53 Côte d’Ivoire
54 Democratic Republic of the Congo
55 Denmark
56 Djibouti
57 Dominica
58 Dominican Republic
59 East Timor
60 Ecuador
61 Egypt
62 El Salvador
63 Equatorial Guinea
64 Eritrea
65 Estonia
66 Eswatini
67 Ethiopia
68 Falkland Islands
69 Faroe Islands
70 Fiji
71 Finland
72 France
73 French Guiana
74 French Polynesia
75 French Southern Territories
76 Gabon
77 Gambia
78 Georgia
79 Germany
80 Ghana
81 Greece
82 Greenland
83 Grenada
84 Guadeloupe
85 Guam
86 Guatemala
87 Guinea
88 Guinea-Bissau
89 Guyana
90 Haiti
91 Heard Island and McDonald Islands
92 Honduras
93 Hong Kong
94 Hungary
95 Iceland
96 India
97 Indonesia
98 Iran
99 Iraq
100 Ireland
101 Isle of Man
102 Israel
103 Italy
104 Jamaica
105 Japan
106 Jersey
107 Jordan
108 Kazakhstan
109 Kenya
110 Kiribati
111 Kosovo
112 Kuwait
113 Kyrgyzstan
114 Laos
115 Latvia
116 Lebanon
117 Lesotho
118 Liberia
119 Libya
120 Liechtenstein
121 Lithuania
122 Luxembourg
123 Macao
124 Macedonia
125 Madagascar
126 Malawi
127 Malaysia
128 Mali
129 Malta
130 Martinique
131 Mauritania
132 Mauritius
133 Mayotte
134 Mexico
135 Micronesia
136 Moldova
137 Mongolia
138 Montenegro
139 Montserrat
140 Morocco
141 Mozambique
142 Myanmar
143 Namibia
144 Nepal
145 Netherlands
146 New Caledonia
147 New Zealand
148 Nicaragua
149 Niger
150 Nigeria
151 Niue
152 North Korea
153 Northern Cyprus
154 Northern Mariana Islands
155 Norway
156 Oman
157 Pakistan
158 Palau
159 Palestine
160 Panama
161 Papua New Guinea
162 Paraguay
163 Peru
164 Philippines
165 Poland
166 Portugal
167 Puerto Rico
168 Qatar
169 Romania
170 Russia
171 Rwanda
172 Réunion
173 Saint Helena
174 Saint Kitts and Nevis
175 Saint Lucia
176 Saint Pierre and Miquelon
177 Saint Vincent and the Grenadines
178 Saint-Martin
179 Samoa
180 Saudi Arabia
181 Senegal
182 Serbia
183 Seychelles
184 Sierra Leone
185 Singapore
186 Sint Maarten
187 Slovakia
188 Slovenia
189 Solomon Islands
190 Somalia
191 South Africa
192 South Georgia and the South Sandwich Islands
193 South Korea
194 South Sudan
195 Spain
196 Sri Lanka
197 Sudan
198 Suriname
199 Svalbard and Jan Mayen
200 Sweden
201 Switzerland
202 Syria
203 São Tomé and Príncipe
204 Taiwan
205 Tajikistan
206 Tanzania
207 Thailand
208 Togo
209 Tonga
210 Trinidad and Tobago
211 Tunisia
212 Turkey
213 Turkmenistan
214 Turks and Caicos Islands
215 Uganda
216 Ukraine
217 United Arab Emirates
218 United Kingdom
219 United States
220 Uruguay
221 Uzbekistan
222 Vanuatu
223 Venezuela
224 Vietnam
225 Virgin Islands, U.S.
226 Western Sahara
227 Yemen
228 Zambia
229 Zimbabwe
230 Åland


Lo que vamos hacer es copiar el objeto paises_r y llamarlo Colombia, luego vamos a reclasificar todos los valores diferentes a 43 con NA (sin valores), es decir escoger Colombia.

Colombia <- paises_r
Colombia[Colombia != 43] <- NA
Colombia
## class       : SpatRaster 
## dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## categories  : NAME_0 
## name        :   NAME_0 
## min value   : Colombia 
## max value   : Colombia
# Compare el código de arriba con este de abajo. ¿Cuales son las diferencias? ¿Podría ser de utilidad?

# Colombia <- paises_r == 43

Visualicemos el mapa de Colombia. Hemos seleccionado lo que queremos, pero el extent: sigue siendo el mismo, es decir todo el globo entero, por eso vemos el país tan pequeño.

plot(Colombia)

Podemos utilizar la función trim() para remover las filas y columnas con valores NA (alrededor de los valores del mapa de Colombia), y como consecuencia se cambia el extent: del nuevo objeto.

Colombia_t <- terra::trim(Colombia)
Colombia_t
## class       : SpatRaster 
## dimensions  : 200, 146, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -79, -66.83333, -4.25, 12.41667  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## categories  : NAME_0 
## name        :   NAME_0 
## min value   : Colombia 
## max value   : Colombia
plot(Colombia_t)

Ahora sí utilicemos la función crop() para extraer los datos de temperatura media en el extent: del objeto Colombia_t.

Colombia_huella <- terra::crop(huella, Colombia_t)
Colombia_huella
## class       : SpatRaster 
## dimensions  : 200, 146, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -79, -66.83333, -4.25, 12.41667  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : wildareas-v3-2009-human-footprint 
## min value   :                              0.00 
## max value   :                             37.86

Visualicemos el nuevo objeto. Como podemos ver la función crop() extrae lo valores de la huella humana con base en el extent: del objeto Colombia_t.

plot(Colombia_huella)

Para tener solamente los valores de la huella humana dentro de los limites de Colombia_t debemos utilizar la función mask(). Veamos:

Colombia_huella_mask <- terra::mask(Colombia_huella, Colombia_t)
plot(Colombia_huella_mask)

Ahora bien, por qué no utilizar la opción mask() desde un principio? Hay varias razones.

  • las dos rásters deben ser del mismo extent:. Al intentar implementar la función mask() directamente al mapa global huella obtenemos el Error. Es decir el mapa huella (es global) y el mapa Colombia_t tienen diferente extent:.
Colombia_huella_2 <- terra::mask(huella, Colombia_t)
## Error: [mask] extents do not match
  • la segunda razón es que es más eficiente para R implementar la función crop(), porque selecciona filas y columnas del área definida por el extent: que queremos. Al contrario, mask() require hacer cálculos más complejos para seleccionar los pixeles que están bordeando la figura del mapa de Colombia.

Tips!

Por eso, es importante implementar primero crop() y luego mask() para acelerar el proceso.

7.6.4. Seleccionando valores en un ráster

En ocasiones necesitamos utilizar solo unos valores determinados dentro del ráster.

Por ejemplo, queremos obtener un mapa con las temperaturas por encima de 25 grados. Para lo cual podemos reclasificar los valores menores a 25 con NA (no values).

tmean_caliente <- tmean_mean
tmean_caliente[tmean_caliente < 25] <- NA
tmean_caliente
## class       : SpatRaster 
## dimensions  : 1080, 2160, 1  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :    mean 
## min value   : 25.0000 
## max value   : 30.9879
plot(tmean_caliente)

Pero tambien podemos necesitar un rango , por ejemplo entre 15 y 20 grados. Para eso utilizamos el operador | que indica “o” (los valores menores a 15 “o” valores mayores a 20, serán iguales a NA).

tmean_templado <- tmean_mean
tmean_templado[tmean_templado < 15 | tmean_templado >20] <- NA
plot(tmean_templado)

Tambien podemos extraer lo valores en una pila de rásters. Por ejemplo, seleccionemos los valore mayores a 10 grados en todas las rásters.

tmean_frio <- tmean
tmean_frio[tmean_frio > 10] <- NA
tmean_frio
## class       : SpatRaster 
## dimensions  : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## names       : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ... 
## min values  :     -45.884,       -44.8,   -57.92575,    -64.1925,    -64.8115,   -64.35825, ... 
## max values  :      10.000,        10.0,    10.00000,     10.0000,     10.0000,    10.00000, ...

7.6.5. Seleccionando valores en un ráster con otro ráster

Vamos a extraer los valores de precipitación global en areas con temperaturas mayores a 20°C (tmean_caliente).

Primero, descarguemos la precipitación utilizando el paquete geogata:

prec <- geodata::worldclim_global(
          var="prec",
          res=10,
          path = "C:/Talleres_R/IntroR_espacial"
                                 )
prec
## class       : SpatRaster 
## dimensions  : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : wc2.1_10m_prec_01.tif  
##               wc2.1_10m_prec_02.tif  
##               wc2.1_10m_prec_03.tif  
##               ... and 9 more source(s)
## names       : wc2.1~ec_01, wc2.1~ec_02, wc2.1~ec_03, wc2.1~ec_04, wc2.1~ec_05, wc2.1~ec_06, ... 
## min values  :           0,           0,           0,           0,           0,           0, ... 
## max values  :         908,         793,         720,        1004,        2068,        2210, ...

Ahora calculemos la precipitación media anual:

prec_mean <-terra::app(prec, mean)
prec_mean
## class       : SpatRaster 
## dimensions  : 1080, 2160, 1  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     mean 
## min value   :   0.0000 
## max value   : 932.5833
terra::plot(prec_mean)

Finalmente, utilicemos la función mask() para extraer los valores de precipitación en areas con temperaturas mayores a 20°C.

prec_caliente <- terra::mask(prec_mean, tmean_caliente)
prec_caliente
## class       : SpatRaster 
## dimensions  : 1080, 2160, 1  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     mean 
## min value   :   0.0000 
## max value   : 932.5833

Visualizemos los dos rásters:

# Creemos dos columnas para los plots
par(mfrow=c(1,2))

# Temperature > 20 °C
plot(tmean_mean*0, legend=F) # Fondo gris para mostrar los continentes 
t_breaks <- seq(minmax(tmean_caliente)[1], minmax(tmean_caliente)[2], 1) # escala de la temperatura
terra::plot(tmean_caliente, add=T, type="interval", plg=list(x="bottom", cex=1.2, title="°C"), breaks = t_breaks)
title("Temperatura media > 20 °C", line = -2, cex.main=1.6)

# Precipitación en áreas con temp  > 20 °C
p_breaks <- seq(minmax(prec_caliente)[1], minmax(prec_caliente)[2], 150)
plot(tmean_mean*0, legend=F)
terra::plot(prec_caliente,  add=T, type="interval", plg=list(x="bottom", cex=1.2, title="mm"), col=rev(topo.colors(15)), breaks=p_breaks)
title("Precipitación en áreas con temp  > 20 °C", line = -2, cex.main=1.6)

7.6.6. Extraer los valores de un ráster

En algún momento de sus análisis necesitará los valores de los rásters para hacer algunos cálculos adicionales o preparar gráficas. El paquete terra ofrece la función values() para extraer lo valores del ráster.

Veamos un ejemplo con la temperatura media y la precipitación global.

Visualicemos una gráfica que muestre los valores pareados de temperatura media y precipitación.

Interesante ver que las precipitaciones medias anuales en áreas frías son menores que áreas más calientes.

plot(values(tmean_mean), values(prec_mean))

Ahora tomemos solamente los datos de precipitación para Colombia.

Colombia_t_prec <- terra::crop(prec_mean,  # precipitación media
                        terra::aggregate(Colombia_t, 2) # agreguemos colombia 2 veces para que sea la misma resolución de la precipitación.
                        )%>%
                  resample(aggregate(Colombia_t, 2))%>% # asegurarnos que las rásters sean iguales en extent, resolution, etc.
                  mask(aggregate(Colombia_t, 2))       # Seleccionar solamente el borde administrativo de Colombia

plot(Colombia_t_prec)

El Pacífico siempre lluvioso!

Seleccionemos también la temperatura media para Colombia.

Colombia_t_tmean <- terra::crop(tmean_mean, aggregate(Colombia_t, 2))%>%
                    resample(aggregate(Colombia_t, 2))%>%
                    mask(aggregate(Colombia_t, 2))
plot(Colombia_t_tmean)

Ahora visualicemos los valores pareados de temperatura media y precipitación para Colombia con relación al global.

Definitivamente tropical!

  plot(values(tmean_mean), values(prec_mean), xlim=c(-55, 30), ylim=c(0, 1000),
       ylab="",yaxt="n", xlab="",xaxt="n", pch=21, col="black")

par(new=T)
plot(values(Colombia_t_tmean), values(Colombia_t_prec), col="red", xlim=c(-55, 30), ylim=c(0, 1000), pch=22)

legend("topleft", legend=c("Colombia", "Global"), col=c("red", "black"), pch=c(22,21), cex=0.8)

7.6.7. Proyectando rásters

Hasta ahora hemos utilizado coordenadas geográficas para estos análisis. Sin embargo, los análisis más locales o regionales requieren de ráster proyectadas. Es decir los datos proyectados en una superficie plana. En especial si queremos hacer cálculos de áreas para informar sobre aspectos del área de estudio (ej., área de bosques fragmentados, etc).

Hay diferentes proyecciones cartográficas, pero la que más vamos a utilizar en el análisis espacial es la que mantienen el sentido del área. Cada país y región tiene proyecciones apropiadas que preservan el area, entonces hay que hacer el ejercicio de buscar la mejor proyección para el área de estudio.

Por fortuna hay muchos sitios para buscar y spatialrefence y epsg.io.

Por ejemplo, si buscamos “Colombia” en spatialreference nos aparecerá algo asi:

En nuestros códigos podemos utilizar cualquiera de los dos formatos ESPG (es el que aparece en el primer pantallazo de ‘spatialreference’) o Proj4.

Por ejemplo si utilizamos EPSG:3116 el equivalente Proj4 es: +proj=tmerc +lat_0=4.5962004166 +lon_0=-74.077507916 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Veamos cómo utilizamos esto en R.

  • Carguemos la cobertura vegetal (solamente la clase árboles, “trees”),
  • seleccionemos los datos para Colombia y
  • utilicemos la proyección EPSG:3116 para reproyectar el ráster.

La función project() de terra incluye varios argumentos:

  • y = hace referencia a “coordinate reference system, (por sus siglas en inglés), en este caso EPSG:3116.
  • method = “bilinear” es util para variables continuas y “near” para variables categóricas. Vamos a utilizar “bilinear porque es la proporción de bosque es cada celda (es decir que es una variable continua).
  • res = hace referencia al tamaño de celda (en metros) que queremos al final. Es decir 1Km por 1Km. Celdas cuadradas, lo que facilitará y mejorará la exactitud (accuracy, en inglés) de los cálculos.
arboles <- geodata::landcover(var= 'trees',
                    path = "C:/Talleres_R/IntroR_espacial"                            )%>%
  terra::crop(disagg(Colombia_t,2))%>%
  terra::resample(terra::disagg(Colombia_t,2))%>%
  terra::mask(disagg(Colombia_t,2))%>%
  terra::project(y="EPSG:3116", method="bilinear", res = 1000)
  
arboles
## class       : SpatRaster 
## dimensions  : 1857, 1358, 1  (nrow, ncol, nlyr)
## resolution  : 1000, 1000  (x, y)
## extent      : 451351.3, 1809351, 19203.33, 1876203  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Colombia Bogota zone (EPSG:3116) 
## source(s)   : memory
## name        : trees 
## min value   :     0 
## max value   :     1
plot(arboles)

Utilicemos la función cellSize() de terra para calcular el área de cada celda de este ráster. Casi todas las celdas son cercanas a 1Km2 (~1,000m x ~1,000m), pero hay algunas muy pequeñas variaciones, en especial en el oriente de Colombia, porque esta proyección esta centrada más en la mitad de Colombia.

Continuemos con nuestro análisis, asumiendo que las celdas son de igual tamaño. Sin embargo, tenga en cuenta reportar la exactitud en sus resultados.

arboles_cellsize <- cellSize(arboles, unit="km")
plot(arboles_cellsize)

7.6.8. Reclasificando rásters

La reclasificación de un ráster es de mucha utilidad en el análisis espacial. Por ejemplo si usted quiere reclasificar la cobertura de árboles en Colombia en cuatro categorias.

Para esto debemos crear una matriz que contenga el valor inicial, el valor final del ráster y que valor queremos al final. Repetir el proceso hasta terminar el rango.

Por ejemplo, definamos tres nuevas categorías de acuerdo a la proporción de árboles y reclasifiquemolas con 1, 2 y 3. En la primera linea se dice, de 0 a 0.25 reclasifiquelo a 1.

m <- c(0, 0.25, 1,
       0.25, 0.5, 2,
       0.5, 1, 3)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
rclmat
##      [,1] [,2] [,3]
## [1,] 0.00 0.25    1
## [2,] 0.25 0.50    2
## [3,] 0.50 1.00    3

Ahora apliquemos esta matriz al ráster.

rc <- classify(arboles,
               rclmat)
plot(rc, legend = FALSE,
     col = c("dark green",  "red", "yellow"), axes = FALSE)
legend("topright",
       legend = c("Alta", "Media",  "Baja"),
       fill = c("dark green", "yellow",  "red"),
       border = FALSE,
       bty = "n") # turn off legend border

Podríamos tener unas datos del área de cada categoría de la siguiente forma:

terra::freq(rc)
##   layer value  count
## 1     1     0     91
## 2     1     1 213491
## 3     1     2 164580
## 4     1     3 761372

7.7. Tiempo de procesamiento

El tiempo de procesamiento cuando se manejan o hacen análísis con rásters incrementa considerablemente por varias razones. Mencionemos algunas:

  • La resolución espacial (tamaño de la celda o grid cell en inglés). Celdas más pequeñas generan mayor tiempo de procesamiento.
  • La cobertura (tamaño de la grilla o extent en inglés). Areas más grandes requieren de mayor tiempo de procesamiento.
  • La complejidad y/o tipo de la función o el análisis que se quiere hacer. Por ejemplo, reporjectar un raster puede tomar mucho tiempo.
  • La combinación de todas las anteriores puede generar mucho tiempo de procesamiento o mucho consumo de memoria.

Veamos un ejemplo para visulalizar cómo se afecta el tiempo de procesamiento cuando manejamos o hacemos cálculos con rásters.

Resolución espacial

Creemos 3 grillas de la misma extensión pero diferente tamaño de celda:

pixel <- c(1, 10, 100, 1000, 5000, 7500, 10000) # Tamaño de celda

pixel_size <- list()
for(i in 1:length(pixel)){
  
pixel_size[[i]] <- rast(ncol=pixel[i], nrow=pixel[i])
  values(pixel_size[[i]]) <- 1:ncell(pixel_size[[i]])

}

Ahora apliquemos una función sencilla a cada ráster y calculemos el tiempo que tarda cada una en procesar.

tiempo_proceso <- list()
tiempo <- list()

for(j in 1:length(pixel_size)){
  

  tiempo_proceso[[j]] <- system.time(sqrt(pixel_size[[j]]))[3]
    
}
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          
tiempo <- unlist(tiempo_proceso)
 
df1 <- data.frame(pixel, tiempo)

df1
##   pixel tiempo
## 1     1   0.00
## 2    10   0.00
## 3   100   0.00
## 4  1000   0.01
## 5  5000   0.22
## 6  7500   2.17
## 7 10000   3.92

Visualicemos en una gráfica para que pueda observar la tendencia.

plot(df1$pixel,df1$tiempo, xlab = "Tamaño de celda", ylab = "Tiempo de procesamiento", )

Por ahora tengamos en cuenta que el tiempo de procesamiento es una limitante cuando manejamos rásters. Las soluciones a este problema van desde la optimización de sus códigos hasta la utlización de supercomputadores para acelerar los procesos.

7.8. Tamaño de los rásters creados por R

Cuand se carga un ráster o se aplica un función a un raster, R trata de utilizar la memoria RAM de su computador. Sin embargo, si los valores del el objeto (es decir un ráster) no pude ser guardado en la memoria RAM, entonces R empieza a crear archivos temporales que generalmente son bastante grandes.

Los archivos temporales se crean cuando no se asigan la ruta y nombre del archivo dentro de la función o en funciones en las cuales no existe esta opción de asignar la ruta y el nombre al archivo final.

Antes de ver unos ejemplos, recordemos que al principio habíamos configurado un folder para guardar los archivos temporales.

## [1] "C:/Talleres_R/IntroR_espacial/Temp_para_borrar/r_tmp_2024-07-21_195156.559877_28580_30265.grd"

Ahora si veamos unos ejemplos. Carguemos un raster (elevación) utilizando el paquete geodata::. La función elevation_global() le permite indicar la ruta

if(!dir.exists("C:/Talleres_R/IntroR_espacial/Input_data")){
dir.create("C:/Talleres_R/IntroR_espacial/Input_data")
}else{
  print("Folder Input_data ya existe")
}
## [1] "Folder Input_data ya existe"
elevation <- geodata::elevation_global(res = 2.5,
                                       path = "./Input_data")
elevation
## class       : SpatRaster 
## dimensions  : 4320, 8640, 1  (nrow, ncol, nlyr)
## resolution  : 0.04166667, 0.04166667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : wc2.1_2.5m_elev.tif 
## name        : wc2.1_2.5m_elev 
## min value   :            -415 
## max value   :            7412
terra::plot(elevation)

Ahora intentemos desagregar este ráster, es decir

elevation_ag <-  terra::aggregate(elevation, fact=10, "mean")

terra::size(elevation_ag)
## [1] 373248
elevation_10 <-  sqrt(elevation)
  

terra::size(elevation_10)
## [1] 37324800
showTmpFiles()
## --- none ---
elevation_agg2 <-  terra::disagg(elevation, fact=2, method="near")
## 
|---------|---------|---------|---------|
===========
                                          
terra::size(elevation_agg2)
## [1] 149299200
showTmpFiles()
## --- none ---
sources(elevation_agg2)
## [1] "C:/Users/jzulo/AppData/Local/Temp/RtmpMVZkU4/spat_UCzrRLqRud378rP_28580.tif"

8. Vectorial